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Abstract — Time delay estimation arises in many applications 
in which a multipath medium has to be identified from pufses 
transmitted through the channef. Various approaches have been 
proposed in the literature to identify time delays introduced by 
multipath environments. However, these methods either operate 
on the analog received signal, or require sampling at the Nyquist 
rate of the transmitted pulse. In this paper, our goal is to develop 
a unified approach to time delay estimation from low rate samples 
of the output of a multipath channel. Our methods result in 
perfect recovery of the multipath delays from samples of the 
channel output at the lowest possible rate, even in the presence 
of overlapping transmitted pulses. This rate depends only on the 
number of multipath components and the transmission rate, but 
not on the bandwidth of the probing signal. In addition, our 
development allows for a variety of different sampling methods. 
By properly manipulating the low-rate samples, we show that 
the time delays can be recovered using the well-known ESPRIT 
algorithm. Combining results from sampling theory with those 
obtained in the context of direction of arrival estimation methods, 
we develop sufficient conditions on the transmitted pulse and the 
sampling functions in order to ensure perfect recovery of the 
channel parameters at the minimal possible rate. Our results 
can be viewed in a broader context, as a sampling theorem for 
analog signals defined over an infinite union of subspaces. 



I. Introduction 

Time delay estimation is an important signal processing 
problem, arising in various applications such as radar [1], 
underwater acoustics [2], wireless communications [3], and 
more. In a typical scenario, pulses with a priori known shape 
are transmitted through a multipath medium, which consists 
of several propagation paths. As a result the received signal is 
composed of delayed and weighted replicas of the transmitted 
pulses. In order to identify the medium, the time delay and gain 
coefficient of each multipath component has to be estimated 
from the received signal. 

In this paper we consider recovery of the parameters defin- 
ing such a multipath medium from samples of the channel 
output. Specifically, we assume that pulses with known shape 
are transmitted at a constant rate through the medium, and 
our aim is to recover the time delays and time-varying gain 
coefficients of each multipath component, from samples of 
the received signal. We focus on the sampling stage, and 
derive methods that ensure perfect recovery from samples 
of the channel output at the minimal possible rate. The 
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proposed sampling schemes are flexible, so that a variety 
of different sampling techniques can be accommodated. Our 
main contribution is the development of efficient sampling 
schemes for the received signal. The resulting sampling rate 
is generally much lower than the traditional Nyquist rate 
of the transmitted pulses, and depends only on the number 
of multipath components and the transmission rate, but not 
on the bandwidth of the transmitted pulse. This can lead to 
significant sampling rate reduction, comparing to the Nyquist 
rate, in applications where only small number of propagation 
paths exists, or when the bandwidth of the transmitted pulse 
is relatively high. This reduction is desirable for practical 
implementation. Sampling at lower rates allows for analog-to- 
digital converters (ADCs) that are more precise (i.e. use more 
bits), and with lower power consumption. In addition, lowering 
the sampling rate can reduce the load on both hardware and 
further digital processing units. 

A classical solution for the time delay estimation problem 
is based on correlation between the received signal and the 
transmitted pulse [1]. However, the time resolution of this 
method is limited by the inverse of the transmitted pulse 
bandwidth. Therefore, this technique is effective only when 
the multipath components are well separated in time, or when 
only one component is present. This approach was originally 
motivated in the analog domain, where the entire analog 
correlation is computed. Performing the correlation in the 
digital domain requires samples of the data at a high sampling 
rate, in order to approximate the continuous correlation. 

In order to resolve closely spaced multipath components, 
various superresolution estimation algorithms were proposed. 
In [4], [5], the MUSIC [6] method was applied in the time 
domain in order to estimate the time delay of each multipath 
component. Hou and Wu [7] were the first to convert the 
time estimation problem to model-based sinusoidal parameter 
estimation, and used an autoregressive method in order to 
estimate the model's parameters. Other works, such as [5], 
[8], [9], relied on the same principle, but different estimation 
algorithms were used: Tufts-Kumaresan SVD algorithm [10], 
TLS-ESPRIT method [11] and a modification of MUSIC, 
respectively. 

In the above superresolution approaches, the sampling stage 
was typically not directly addressed. Most of these works rely 
on uniform pointwise samples of the received signal, at a high 
sampling rate. In [7] and [9] the required sampling rate is the 
Nyquist rate of the transmitted pulse. Since often the pulses are 
chosen to have small time support, the bandwidth can be quite 
large, corresponding to a high Nyquist rate. In [8] and in the 
frequency domain algorithm proposed in [5], the sampling of 
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the received analog signal is not mentioned explicitly. Since 
these algorithms involve operations in the analog frequency 
domain, effectively they also require sampling at the Nyquist 
rate. The time domain algorithms proposed in [4] and [5] can 
theoretically recover the time delays from a low sampling rate, 
which depends on the number of multipath components. How- 
ever, the sampling considerations were not directly addressed 
in these works, and no concrete conditions on the transmitted 
pulse and the samples were given, in order to ensure unique 
recovery of the delays from the samples. 

Besides the sampling stage which is not studied in previous 
works, another assumption underlying all the methods above is 
that the receiver has access to a single experiment ([7], [8], [9]) 
or multiple non-overlapping experiments ([4], [5], [9]) on the 
channel. In each experiment a pulse is transmitted through the 
multipath medium, and it is required that all the returns vanish 
before the next experiment is conducted. This imposes the 
constraint that the transmitted pulse is sufficiently time limited, 
which can be problematic in certain scenarios. For example, 
in wireless communications, modulated pulses are transmitted 
at a constant symbol rate through the medium. In this case, we 
cannot consider the observed signal over one symbol period as 
an independent experiment, since it will generally be affected 
by reflections caused by adjacent symbols. 

In Section [TT] we propose a general signal model, that can 
describe the received signal from a time-varying multipath 
medium. An advantage of our model is that it does not require 
the assumption of non-overlapping experiments, and allows 
for general pulse shapes. We then formulate the medium 
identification problem as a sampling problem, in which the set 
of parameters defining the medium have to be recovered from 
samples of the received signal at the lowest possible rate. To 
this end we develop a general sampling scheme, which consists 
of filtering the received signal with a bank of p sampling 
filters and uniformly sampling their outputs. This class of 
sampling methods is common in sampling theory [12] and can 
accommodate a wide variety of sampling techniques, including 
pointwise uniform sampling. Given K multipath components, 
we show that at least 2K sampling filters are required in order 
to perfectly recover the time delays. We then develop explicit 
sampling strategies that achieve this minimal rate. In particular, 
we derive sufficient conditions on the transmitted pulse and the 
choice of sampling filters, which guarantee unique recovery of 
the channel parameters. 

In order to recover the channel parameters from the given 
samples we combine results from standard sampling theory, 
with those of direction of arrival (DOA) algorithms [6], [11], 
[13], [14]. Specifically, by appropriate manipulation of the 
sampling sequences, we show that we can formulate our 
problem within the framework of DOA methods. We then rely 
on the estimation of signal parameters via rotational invariance 
techniques (ESPRIT) [11], developed in that context. The 
unknown delays are recovered from the samples by first 
applying a digital correction filter bank, and then applying the 
ESPRIT algorithm. Once the time delays are identified, the 
gain coefficients are recovered using standard sampling tools. 

The sampling schemes we develop for the channel identi- 
fication problem treated in this paper, can be viewed in the 



broader context of sampling theory. In Section [V] we discuss 
in detail the relationship between our problem and previous 
related setups treated in the sampling literature: sampling 
signals from a union of subspaces [15], [16], compressed 
sensing of analog signals [17], [18], [19], [20], [21], [22], 
and finite-rate of innovation (FRI) sampling [23], [24]. The 
results we develop here can be viewed as a special case of 
analog compressed sensing over an infinite union of spaces, 
and therefore extends previous work, which focused on finite 
unions, to a broader setting. In comparison with FRI tech- 
niques, our approach allows for lower sampling rates and at 
a lower computational cost. Furthermore, we do not need to 
impose stringent conditions on the pulse shapes, as required 
by FRI methods. 

This paper is organized as follows. In Section [TT1 we 
describe our signal model. A general sampling scheme of 
the received signal is proposed in Section [III] Section [IV] 
describes the recovery of the unknown delays from the sam- 
ples, and provides sufficient conditions ensuring a unique 
recovery. Relation to previous work is discussed in detail in 
Section [V] Section [VI] describes an application example of 
channel identification in wireless communications. Numerical 
experiments are described in Section IVHI 

II. Signal Model and Problem formulation 

A. Notations and Definitions 

Matrices and vectors are denoted by bold font, with lower- 
case letters corresponding to vectors and uppercase letters to 
matrices. The ?ith element of a vector a is written as a„, and 
the ijth element of a matrix A is denoted by Ay . Superscripts 
(•)*, (-) T and (-) H represent complex conjugation, transpo- 
sition and conjugate transposition, respectively. The Moore- 
Penrose pseudo-inverse of a matrix A is written as A^. The 
identity matrix of size n is denoted by I„. 

The Fourier transform of a continuous-time signal x (t) £ 
Li is defined by X [w) — J™ x (t) e~ JW *dw, and 

/oo 
x(t)y*(t)dt, (1) 
-oo 

denotes the inner product between two continuous-time sig- 
nals. The discrete-time Fourier transform (DTFT) of a se- 
quence a [n] £ £2 is defined by 

A(e juT ) =5>[n]e-^" T , (2) 
and is periodic with period 2ir/T. 

B. Signal model 

We consider the class of signals that can be written in the 
form 

K 

x (t) = ^2 ^2 a k M 9 (t-tk- nT) , (3) 

fc=l riGZ 

where t = {ifc}k =1 is a set of K distinct unknown time delays 
in the continuous interval [0, T), ak [n] is an arbitrary sequence 
in £2, and g (t) £ L2 is a known function. Each signal from 
this class can describe the propagation of a pulse with known 
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shape g (t) which is transmitted at a constant rate of 1/T 
through a medium consisting of K paths. Each path has a 
constant delay tk, and a time-varying gain, which is described 
by the sequence a k [n]. In cases where the transmitted pulses 
are amplitude modulated, the sequence dk [n] describes the 
multiplication between the pulse amplitude and the gain co- 
efficient of each path. In Section [VI] we will discuss more 
thoroughly an example of a communication signal transmitted 
through a multipath time-varying channel. 

Our problem is to determine the delays r and the gains 
ctk[n] from samples of the received signal x(t), at the minimal 
possible rate. Since these parameters uniquely define x(t), our 
channel identification problem is equivalent to developing ef- 
ficient sampling schemes for signals of the form ||3}, allowing 
perfect reconstruction of the signal from its samples. 

The model (f3]) is more general than that described in 
previous work [7], [4], [5], [8], [9], which is based on single or 
multiple experiments on the medium. In each experiment the 
received signal is observed over a finite time window, which 
is synchronized to the transmission time of the pulse. More 
precisely, the received signal in the nth experiment is given 
by 

K 

x n (t) = J2 a kn9(t-t k ), t£[0,T), (4) 

k=l 

where tk is the delay of the fcth multipath component which is 
constant in all the experiments, and ctkn is the gain coefficient 
of the kth multipath component at the nth experiment, which is 
generally varying from one experiment to another. This model 
can be seen as a special case of ([3]) with additional constraints 
on the pulse g (t) and the transmission rate 1 /T. Indeed, we 
can write the received signal on the nth experiment as 

x n (t) =x(t- nT) t e [0, T) , (5) 

if we require that 

g (t - t k - nT) = for all n ^ 0. (6) 

This requirement means that the pulse g (t) has finite time 
support, and that the repetition period of the pulses T, is long 
enough such that all the reflections from one pulse vanish 
before the next pulse is transmitted. On the other hand, our 
signal model does not require these constraints, so that it can 
support infinite length pulses and allows interference between 
experiments. 

III. Sampling Scheme 

A. Known delays 

Before we treat our sampling problem of signals of the 
form (O, we first consider a simpler setting in which the 
delays tk are known. In this case our signal model is a special 
case of the more general class of signals that lie in a shift- 
invariant (SI) subspace. For such signals classes, there are well 
known sampling schemes that guarantee perfect recovery at the 
minimal possible rate [12], [17], [25]. Below, we review the 
main results in this setting which will serve as a basis for our 
development in the case in which the delays are unknown. 



A finitely-generated SI subspace of is defined as [26], 
[27], [28] 

a = < x (t) = ^2 ^2 ak ["] 9k (* ~ nT ) : ° fe n e ^ > ■ 

t k=ln£Z ) 

(7) 

The functions g). (t) are referred to as the generators of A. 
In order to guarantee a unique stable representation of any 
signal x(t) <E A by coefficients a k [n], the generators g k (t) 
are typically chosen to form a Riesz basis of A [26], [27]. 
Clearly, the signal model in (O is a special case of (jT) with 
generators gk(t), obtained from K delayed versions of g (t): 

g k (t)=g{t-tk), l<k<K. (8) 

One way to sample a signal of the form (O is to use K 
parallel sampling channels [17]. In each channel the signal 
is first filtered with an impulse response s» (— t), and then 
sampled uniformly at times t = nT to produce the sampling 
sequence ci [n], as depicted in the left-hand side of Fig. 03 
The sampling sequence at the output of the Ah channel can 
be written as 

c, [n] = (x (t) ,s t (t- nT)) , 1 < I < K. (9) 

By analyzing the DTFTs Ce (e :,wT ) of the sequences 
C*[ti],1 < £ < K, it was shown in [17] that the sequences 
o-t [n] , 1 < I < K, which define the signal x (<), can be 
recovered from ct[n] using an adequate multichannel filter 
bank. Specifically, let c (e 3ulT ) , a {e ju)T ) denote the length-if 
column vectors whose Mi elements are Ci (e JwT ) , Ai (e jajT ) 
respectively. Then, it can be shown that 

c (e^ T ) = M SG (e^ T ) a (e^ T ) , (10) 

where M.$g (e^ 7 ) is a K x K matrix with £kx\\ element 

*s eGk (e^ T ) = ^Y,StL- pn) G k L - . 

(ID 

Here Gk (w) and Si (uS) denote the Fourier transforms of gk (t) 
and S£ (t) respectively. If this matrix is stably invertible a.e in 
u>, then the sequences an [n] can be recovered from the samples 
using the multichannel filter bank whose frequency response 
is given by M^ G (e JwT ), as depicted in the right-hand side 

of Fig. m 

The proposed sampling scheme achieves an average sam- 
pling rate of K/T since there are K sampling sequences each 
at rate 1/T. Intuitively, this approach requires one sample per 
degree of freedom of the signal x (t): for a signal of the form 
(O, under the assumption that the time delays are known, or 
a signal of the form ((7J, in each time period T there are K 
new parameters. 

B. Unknown delays 

We now address our original problem of designing a sam- 
pling scheme for signals of the form ([3]) with unknown delays. 
We propose a system similar to that used in the case of 
known delays, comprised of parallel sampling channels. Since 
now there are more degrees of freedom in the signal x (t), 
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a! (-() ?\ C 'W 

t = nT 



4H 



Fig. 1. Sampling and reconstruction scheme for the case of known delays 



intuitively we will require at least the same number of channels 
as in the case of known delays. Denoting the number of 
channels by p, this implies that p > K. As we will see, 
under certain conditions on the sampling filters and pulse g(t), 
p = 2K sampling filters are sufficient to guarantee perfect 
recovery of x(t) from the samples. We will also show that this 
is the minimal possible rate achievable for all signals x(t). 

In each channel of our sampling system the signal x (t) 
is pre-filtered using the filter s| (— t) and sampled uniformly 
at times t = nT resulting in the samples {9j. In the Fourier 
domain, we can write (0 as 



mez, 



(12) 

From the definition of x(t), its Fourier transform can be 
written as 

K 



k=l nGZ 
K 

= ^2A k {e^ T )G(cj)e-^, (13) 



k=l 



where A k (e J " T ) denotes the DTFT of the sequence a k [n], 
and G (uj) denotes the Fourier transform of g (t). Substituting 
([T3l into (O, we have 

Kirn = EM*n?r?: 



I.— i ^ roGZ ^ ^ 



= E^( e ' wT K^4£^- % 



k=l 



Here M (e^ 7 ', r) is a p x K matrix with £fcth element 



2tt 

T 1 



■G 



2tt 



uj m I e J T mt ' 



(16) 



and D ( 



t) is a diagonal matrix with fcth diagonal 



element equal to e i utk . Defining the vector b (i 

b(e^ T ) =D (e^ T ,r)a(e^ T ), 
we can rewrite (fl3T l in the form 



as 



) = M (e^ T , r) b (e 



(17) 



(18) 



Our problem can then be reformulated as that of recovering 
b (e^ 7 ) and the unknown delay set r from the vectors 
c(e Ja,T ), for all oj e [0,^). Once these are known, the 
vectors a (e : > ulT ) can be recovered using the relation in (T7\ . 

To proceed, we focus our attention on sampling filters Se(u>) 
with finite support in the frequency domain, contained in the 
frequency range 

"2tt 2tt , 

"7, 7F VP + 7) > 



T T 



(19) 



where 7 6 Z is an index which determines the working 
frequency band T . This choice should be such that it matches 
the frequency occupation of g (t) (although git) does not have 
to be bandlimited). This freedom allows our sampling scheme 
to support both complex and real valued signals. Under this 
choice of filters, each element M (e?^ , r) of ( TToT i can be 
expressed as 

p 

M ek (e^ T ,r) = E W, m (e^ T ) N„ ife (r) , (20) 

m— 1 

where W (e- JwT ) is a p x p matrix whose ^mth element is 
given by 



We 



^s* e (u + y ( m - 1 + 7) 



2tt 

■G(w + — (m-1 + 7) 



(21) 



2tt 



■G [ w - — m I e J - mtfc : 



(14) 



where the first equality is a result of the fact that Af. (e- 7 " 7 ) 
is 27r/T-periodic. 

From now on, we will assume that oj e [0, and all the 
expressions in the DTFT domain are 2n/T periodic. Denoting 
by c (e 3 " T ) the length-p column vector whose £th element is 
Ct (e : ' ajT ) and by a (e^ T ) the length-if column vector whose 
fcth element is A k (eJ uT ), we can write ([Pil l in matrix form 
as 

c (e^ T ) = M (e^ T , r) D (e^ T , r) a (e^ T ) . (15) 



and N (r) is a p x K Vandermonde matrix with mfcth element 
N mfe (r) = e- i ^( m - 1 +^. (22) 
Substituting (f20]l into ( fT8l , 

c(e^ T ) = W (e^ T ) N (r) b(e^ T ). (23) 

If W (e J " T ) is stably invertible, then we can define the 
modified measurement vector d (e : * ulT ) as 



d (e^ T ) = W- 1 (e^ T ) c (e^ T ) . 

This vector satisfies 

d(e^ T ) =N(r)b(e^ T ). 

Since N (r) is independent of w, from the linearity of the 
DTFT, we can express (l25T l in the time domain as 

d [n] = N (r) b [n] , neZ. (26) 



(24) 
(25) 
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Fig. 2. Sampling and reconstruction scheme for the case of unknown delays 



The elements of the vectors d [n] and b [n] are the discrete 
time sequences, obtained from the inverse DTFT of the 
elements of the vectors b ie^ T ) and d (e- 7 " 2 ") respectively. 

Equation (|25T > and its equivalent time domain representation 
(|26"1 >. describe an infinite set of measurement vectors, each ob- 
tained by the same measurement matrix N (t), which depends 
on the unknown delays r. This problem is reminiscent of the 
type of problems that arise in the field of DOA estimation 
[14], as we discuss in the next section. One class of efficient 
methods for DOA recovery, are known as subspace methods 
[14]. These techniques have subsequently been applied to 
many other problems such as spectral estimation [29], system 
identification [30] and more. Our approach is to rely on these 
methods in order to first recover r from the measurements. 
After r is known, the vectors b (e :,ajT ) and a [eP^) can be 
found using linear filtering relations by 



b (e^ T ) = N 1 " (r) d i 



(27) 



Since N (r) is a Vandermonde matrix, its columns are linearly 
independent, and consequently N^N = Ik- Using (fTTI l. 

r) N 1 " (r) d (e^ T ) . (28) 



ale 



D 1 (e^ T , 



The resulting sampling and reconstruction scheme is depicted 
in Fig. |2] 

Our last step, therefore, is to derive conditions on the 
filters s| (— t) and the function g (t) in order that the matrix 
W {e-' u)T ) is stably invertible. To this end, we can decompose 
the matrix W (e jujT ) as 

W (e^ T ) = S (e^ T ) G (e^ T ) (29) 

where S (e^ 7 ) is a p x p matrix with ^mth element 

St m (e^ T ) = is? L + y (™ - 1 + i)) (30) 

and G (e JwT ) is a p x p diagonal matrix whose 771th diagonal 
element is given by 

G mm (e^ T ) = G (w + ^ (ro - 1 + 7)) . (31) 

Each one of these matrices needs to be stably invertible. 
Therefore, from QTT l the condition that the function g (t) needs 
to satisfy is that 



< a < \G (w)| < b < 00 a.e to £ T. 



(32) 



In addition the filters s\ (— t) should be chosen in such a way 
that they form a stably invertible matrix S (e Ja,T ). Examples 
of such filters are given in the next subsection. 



We note here that the conditions given above guarantee 
a stable digital correction filter bank W _1 (e JajT ), however 
generally it will be comprised of infinite length digital filters. 
Practical implementation of these filters can be achieved by 
truncating the impulse response. The length of the resulting 
filters will affect the total delay of our proposed method, which 
will generally be longer than that of the methods described in 
Section HI due to the additional digital filtering stage. 

We summarize the results so far in the following proposi- 
tion. 

Proposition 1: Let ct [n] = (x (t) ,se(t — nT)) , 1 < £ < 
p be a set of p sequences obtained by filtering the signal 
x (t) defined by OJ with p filters s| (— t) and sampling 
their outputs at times nT. Let S^(w) be supported on T = 
[2e 7) 22 (p + 7 )], and let fi = [0, 22). If the function g (t) 
satisfies the condition in ( l32t and the matrix S (e JwT ), defined 
by d30"l >. is stably invertible a.e uj £ Q, then the delays r and 
vector b (e^" 7 ) can be found from the set of equations 



N(r)b 



(33) 



using subspace methods, described in the next section. Here 
N (t) is a p x K Vandermonde matrix with mfcth element 
e -3'^-(m-l+7)** anc [ 



d (e^ T ) = W- 1 ( 



-1 (j<T\ c ( e ^ T 



(34) 



with W (e JwT ) defined by (f2TT >. The sequences ai[n] can then 
be recovered by 

a (e^ T ) = D- 1 (e JwT ,T) N^d (e^ T ) , to £ fi (35) 

where D(e Ja ' T,r ) is a diagonal matrix with diagonal elements 



C. Examples of filters 

We now provide some examples of filters si (t) satisfying 
the required conditions. 

1) Complex bandpass filter-bank: The first example is a set 
of complex bandpass filters. We assume that the working band 
is T = [0, %jtp] (7 = 0), and the function g (t) satisfies (|32l ) 
on that frequency range. We choose the filters s* t (— t) as ideal 
bandpass filters, covering consecutive frequency bands: 



S e (oo) 



to £ [(£-1) 
otherwise. 



2tt /> 271-1 
T ■>*■ T } 



(36) 



The resulting matrix S (e JwT ) is diagonal, and stably invert- 
ible. This example generalizes to any valid working band, 
given by ( fr9l , by shifting the frequency response of the filters. 

We now provide an example demonstrating the importance 
of the sampling filter. 

Example 1: We consider the case where g(t) — 8 it) and 
there are K = 2 diracs per period of T = 1, as illustrated in 
Fig. 0a). The sampling scheme described above, consisting 
of a complex bandpass filter-bank, is used. In Figs. |3jb)- 
(d), we show the outputs of the first 3 sampling channels. 
This example demonstrates the need for the sampling filters 
when sampling short-length pulses at a low sampling rate. 
The sampling kernels have the effect of smoothing the short 
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(a) 



(b) 



pass filter with cutoff irp/T, followed by a uniform sampler 
at a rate of p/T. 

IV. Recovery of the Unknown Delays 

We have seen in the previous section that perfect recon- 
struction of a signal x (t) of the form (0, is equivalent to that 
of recovering the delays t from the modified measurements of 
d2oT ). As we now show, this problem is similar to that arising 
in DOA estimation. 






(O 



(d) 



Fig. 3. Stream of diracs. (a) K = 2 diracs per period T = 1. (b)-(d) 
The outputs of the first three sampling channels. The dashed lines denote the 
sampling instants. 



pulses (diracs in this example). Consequently, even when the 
sampling rate is low, the samples contain information about the 
signal. In contrast, if we were to sample the signal in Fig. a) 
directly at a low rate, then we would often obtain only zero 
samples which contain no information about the signal. 

2) Delayed channels: In this example we assume p is even 
and define the working band as 



(37) 



(7 — — f>/2). We also assume that g(t) satisfies d32| l. We 
choose the £th filter as a delay of Ag £ [0, T) followed by an 
ideal low pass filter. Thus, 




otherwise. 



(38) 



With this choice of filters, the £mth element of the matrix 
S (e jujT ) defined in d30l ) is given by 

= e j( IJ -T(p/ 2 )) A ' e JT( m - 1 ) A ' i (39) 



The matrix S (e^ 7 ) can be expressed as 
S (e^ T ) = * (e^ T ) F, 



(40) 



where <fr ( e,J " T ) lS a diagonal pxp matrix whose £th diagonal 
element is e^" J ~~ < " p ^ 2 ^ Ai , and F is a Vandermonde matrix 
whose £mth element is given by e i^-( m - 1 ) A f From (|40l > it 
can be seen that S (e J " T ) is invertible for all w 6 fi, when 
the delays in each channel Ag are distinct. 

One special case of this choice of sampling filters is when 
the delays are uniformly spaced, i.e Ag — (£ ~ 1) T/p. In this 
case our sampling scheme can be implemented by an ideal low 



A. Relation to direction of arrival estimation 

In DOA estimation [6], [11], [13], [14], K narrow band 
sources impinge on an array, composed of p sensors, from 
distinct DOAs. The goal is to estimate the DOAs of the sources 
from a set of M measurements, obtained from the sensors 
outputs at distinct time instants. 

The DOA estimation problem can be formulated using the 
following measurement model 



X = A(6)S 



(41) 



where X is a p x M matrix, composed of the measurements 
in its columns, S is a K x M matrix consisting of the sources 
signals in its columns and A (9) is a p x K matrix which 
depends on the set of unknown DOAs 8. The structure of 
A (O) is such that its fcth column, denoted a(#fe), depends 
only on the DOA Ok of the fcth source. The vector a(9k) 
is referred to as the steering vector of the array toward 
direction 9 k- The set containing all possible steering vectors, 
i.e, {a (9) , 8 G [0, 2tt)} is referred as the array manifold. 
Given X, the problem is to recover the DOAs 9k, and the 
sources S. 

The set of equations in d26l > has the same form as ((4111 . 
The fcth column of the matrix N (r) depends only on the fcth 
unknown delay tk, and can be described by the vector n (tk), 
where 

n (t) = [ e-^ 7 * e -^( 1+ ' 7 ) t ■•■ e"^ (p_1+7) * ] T . 

(42) 

The array manifold in our setting is the set of vectors 
{n(t),t G [0, T)}. Therefore, we can adopt DOA methods 
to estimate the unknown delays. The only difference between 
the two problems is that in our setting, we have infinitely 
many measurement vectors, in contrast to the DOA problem in 
which X has finitely many columns. This will require several 
adjustments, which we will detail in the ensuing subsections. 

Two prominent methods used for DOA estimation are 
MUSIC (Multiple Signal Classification) [6] and ESPRIT (Es- 
timation of Signal Parameters via Rotational Invariance Tech- 
niques) [11]. These algorithms belong to a class of techniques, 
known as subspace methods, which are based on separating 
the space containing the measurements into two subspaces, the 
signal and noise subspaces. Estimating the unknown set of pa- 
rameters using MUSIC involves a continuous one-dimensional 
search over the parameter range. This procedure can be costly 
from a computational point of view. The ESPRIT approach 
can estimate the unknown set of parameters more efficiently, 
by imposing the additional requirement that the measurement 
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matrix is rotationally invariant. We describe this property in 
subsection II V-CI and show that in our case N (r) satisfies this 
condition, and therefore we use the ESPRIT approach. 

We note that although the MUSIC and ESPRIT meth- 
ods were originally developed as DOA estimators, these ap- 
proaches and other subspace methods have been used success- 
fully in other fields. 

B. Sufficient conditions for perfect recovery 

We now rely on results obtained in the context of DOA 
estimation in order to develop sufficient conditions for a 
unique solution to ( f26b . Such a solution consists of the infinite 
set of vectors b [n] , n e Z and the unknown delays r. 

Conditions for a unique solution (0, S) for (|4TT > where 
derived in [31]. Since [31] deals with a finite number of 
measurements, we have to extend the results to our case, which 
consists of an infinite number of measurements. The analysis 
in [31] requires a preliminary condition that any subset of p 
distinct steering vectors from the array manifold is linearly 
independent. In our case this condition translates into the 
requirement that any set of p vectors n(f,-),l < i < p 
associated with distinct delays U £ [0, T) , 1 < i < p are 
linearly independent. From (|42j, any such set forms a p x p 
Vandermonde matrix, and are therefore linearly independent. 
Therefore, this condition automatically holds in our problem 
without forcing any additional constraints. 

To derive sufficient conditions for a unique solution of the 
set of infinite equations d26l i we introduce some notation. We 
define the measurement set d [A] as the set containing all 
measurement vectors d [A] = {d [n],n£ Z}. Similarly, we 
define the unknown vector set b [A] as b [A] = {b[n],n£ Z}. 
We may then rewrite (l26l l as 

d [A] = N (r) b [A] . (43) 

The following proposition provides sufficient conditions for a 
unique solution to d43l . 

Proposition 2: If (f,b [A]) is a solution to d43l . 

p > 2 J ftT-dim(span(b[A])) (44) 

and 

dim (span (b [A])) > 1, (45) 

then (r, b [A]) is the unique solution of ( [43) . 
The notation span (b [A]) is used for the minimal dimension 
subspace containing the unknown vector set b [A]. The condi- 
tion ( f43T > is needed to avoid the case where b [A] = 0. In this 
case clearly the set r can not recovered uniquely. 

Proof: We denote r = dim (span (b [A])). From (|45T > 
r > 1. Therefore, there exist a finite subset A — {rii}l_-. C A, 
such that the vector set b[A] spans an r-dimensional subspace: 

dim (spun (b[A]^ = r. (46) 

By defining the matrices B and D as the matrices consisting 
of the vector sets b[A] and d[A], we can write 

D = N (f) B. (47) 



Clearly, from its construction, the rank of the matrix B is r. 
From (@D, 

p > 2K - r. (48) 

According to Theorem 1 in [31], the solution (f, B) is the 
unique solution to d47j > under the condition ( 148) . 

Since the set of unknown delays f is the unique solution to 
the finite set of equations ( |47| i, it is also a unique solution to the 
infinite set of equations d43l . Once f is uniquely determined, 
the matrix N (f) is known. Since every vector of the vector 
set b[A] is contained in C K , 

dim (span (b[A]\ ) < K. (49) 

Therefore, according to ( f44l > p > K. The matrix N (f) is 
a p x K Vandermonde matrix which consist of K linearly 
independent vectors. Therefore, for every n 6 A, if b [n] is a 
solution to 

d [n] = N (f) b [n] (50) 

then it is the unique solution. ■ 
Proposition [2] suggests that a unique solution to the set of 
equations d26i i is guaranteed, under proper selection of the 
number of sampling channels p. This parameter, in turn, 
determines the average sampling rate of our sampling scheme, 
which is given by p/T. The condition d44i i depends on the 
value of dim (span (b [A])), which is generally not known in 
advance. According to our assumption dim(span(b[A])) > 1, 
therefore in order to satisfy the uniqueness condition d44i > for 
every signal of the form (0, we must have p > 2K — 1 
sampling channels or a minimal sampling rate of 2K/T. 
Comparing this result to the minimal sampling rate in the case 
when the delays are known in advance, there is a penalty of 
2 in the minimal rate. 

In Section IV-AI we show that our signal model, described 
in (|3), can be considered as part of a more general framework 
of signals that lie in a union of SI subspaces [15]. It was 
shown in [15] that the theoretical minimum sampling rate 
required for perfect recovery of such a signal from its samples 
is 2K/T. Therefore, according to the results of Proposition [2] 
our sampling scheme can achieve the minimal sampling rate 
required for signals of the form (0. 

The minimal sampling rate of 2K/T, which is achieved 
by our scheme, does not depend on the bandwidth of the 
pulse g (t), but only on the number of propagation paths 
K. In applications where the number of propagation paths 
is relatively small, or the bandwidth of the transmitted pulse 
is high, our approach can provide a sampling rate lower 
than the Nyquist rate. More precisely, when 2K/T < W, 
where W is the bandwidth of transmitted pulse, our method 
can reduce the sampling rate relatively to the Nyquist rate. 
For example, the setup in [32], used for characterization of 
ultra-wide band (UWB) wireless indoor channels, consists of 
pulses with bandwidth of W = 1GHz transmitted at a rate 
of 1/T = 2MHz. Under the assumption that there are 32 
significant multipath components, our method can reduce the 
sampling rate down to 128MHz compared with the 2GHz 
Nyquist rate. 
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Beside the theoretical interest, sampling rate reduction is 
also important for implementation considerations. For practical 
ADCs, which perform the sampling process, there is a trade- 
off between sampling rate and resolution [33]. Therefore, 
reducing the sampling rate allows the use of more precise 
ADCs, which can improve the time delay estimation. The 
power consumption of an ADC can also be reduced by 
lowering the sampling rate [33]. In addition, a lower rate leads 
to more efficient digital processing hardware, since a smaller 
number of samples has to be processed. This also allows 
performing the digital processing operations in real time. 

C. Recovering the unknown delays 

According to Proposition [2] in order be able to perfectly 
reconstruct every signal of the form ©, our sampling scheme 
must have p > 2K sampling channels. We assume throughout 
that this condition holds. 

We now describe an algorithm for the recovery of the un- 
known delays from the measurement set d [A], which is based 
on the ESPRIT [11] algorithm. One of the conditions needed 
in order to use the ESPRIT method is that the correlation 
matrix 

R bb = b [n] h H [n] , (51) 

is positive definite. In order to relate this condition to our 
problem, we state the following proposition from [21]. 

Proposition 3: If the sum (l5TT i exists, then every matrix V 
satisfying R bb = YV H has column span equal to span (b [A]). 

An immediate corollary from Proposition [3] is that R^ >~ 
is equivalent to the condition dim (span (b [A])) = K. In this 
case, which we refer to as the uncorrelated case, we can apply 
the ESPRIT algorithm on the measurement set d [n] in order to 
recover the unknown delays. The case dim (span (b [A])) < K, 
will be referred to as the correlated case. In this setting the 
condition R^ y doest not hold, and the ESPRIT algorithm 
cannot applied directly. Instead, we will use an additional stage 
originally proposed in [34], [35]. 

Note, that 

R dd = J]d[n]d H [n] 

= N(r)P>>[n]b ff [n]W(T) 

VnGZ / 

= N (t) K bb N H (r) . (52) 

Since for any set of delays the matrix N (r) has full 
column-rank, the ranks of the matrices Rdd and R&& are equal. 
Therefore, the decision whether we are in the uncorrelated or 
correlated case can made directly from the given measure- 
ments by forming the matrix Rdd- 

1) Uncorrelated Case: From d52l , under the assumption 
that the matrix R^ is positive definite, it can be shown that 
the rank of the matrix Rdd is K. Moreover, the matrices Rdd 
and N (r) have the same column span which is referred as 
the signal subspace. By performing a singular value decom- 
position (SVD) of the matrix Rdd, we can obtain K vectors, 
which span the signal subspace, by taking the K left singular 



vectors associated to the non-zero singular values of Rdd- We 
define the p x K matrix E s as the matrix containing those 
vectors in its columns. 

Now, we will the exploit the special structure of the Vander- 
monde matrix. We denote the matrix Nj (r) as the sub matrix 
extracted from N (r) by deleting its last row. In the same way 
we define Nf (r) as the sub matrix extracted from N (r) by 
deleting its first row. The Vandermonde matrix N (r) satisfies 
the following rotational invariance property: 

N T (r) =N x (t)R(t) (53) 

where R (r) is a diagonal K x K matrix, whose fcth diagonal 
element is given by Rkk (t) = e~ j27rtfc/ ' T . Since the matrices 
N (t) and E s have the same column span, there exists an 
invertible K x K matrix T such that 

N (r) = E S T. (54) 

By deleting the last row in d54l i we get 

N l (r) = E si T. (55) 

Similarly, deleting the first row in (l54l and using the rotational 
invariance property (l53l l. we have 

N i (r)R(r)=E sT T. (56) 

Combining d55l > and (l56l > leads to the following relation 
between the matrices E B f and E s ^: 

E sT =E sJ TR(r)r 1 . (57) 

The matrix E s f is a (p — 1) x K (p > K) matrix with full 
column rank. Therefore E^, E s j = Ik- Using ( TSTT i we define 
the following K x K matrix $ as 

$ = E| 1 E st =TR(r)r 1 . (58) 

From ( 1581 it is clear that the diagonal matrix R (r) can be 
obtained from the matrix $ by performing an eigenvalue 
decomposition. Once the matrix R (r) is known, the unknown 
delays can be retrieved from its diagonal elements as 

i fe = -— arg(R fcfc (r)). (59) 

In summary, our algorithm consist of the following steps: 

1) Construct the correlation matrix Rdd = 
En^d[n]d H [n]. 

2) Perform an SVD decomposition of R^ and construct 
the matrix E s consisting of the K singular vectors asso- 
ciated with the non-zero singular values in its columns. 

3) Compute the matrix $ = E^E s j. 

4) Compute the eigenvalues of <I>, Aj, i = 1, 2, . . . , K. 

5) Retrieve the unknown delays by = — ^arg(Ai). 

2) Correlated Case: When the condition R b b >- is not 
satisfied the ESPRIT algorithm cannot be applied directly on 
the vector set d [A]. In this case the rank of Rdd is smaller 
than K, and therefore its column span is no longer equal 
to the entire signal subspace. To accommodate this setting, 
we perform an additional stage before applying the ESPRIT 
method, based on the spatial smoothing technique proposed in 
[35], [34]. 
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To proceed, we define M = p — K length- (K + 1 ) sub 
vectors 

di [n] = [ di [n] d i+ i[n] ... d i+K [n]] T . (60) 
We define the smoothed correlation matrix Udd as 
- 1 M 

^ EE d *w<*f w. (6i) 

i=l nGZ 

Under our assumptions p > 2K, therefore M > K. 
According to [35], when M > K the rank of the smoothed 
correlation matrix is K regardless of the rank of the matrix 
Hbb- We will refer now to column rank of Rdd as the signal 
subspace, and can then apply the ESPRIT algorithm on this 
matrix. 



provided, these do not necessarily imply that there exists an 
efficient recovery algorithm, which can recover the signal from 
its samples at the minimal rate. Our aim in this work, is 
to provide concrete recovery techniques, that are simple to 
implement, for signals over an infinite union of SI subspaces. 

In summary, in this work we focus on a special case of 
signals that lie in an infinite union of SI subspaces. For this 
case, in contrast to [15], we provide a concrete reconstruction 
method. This method achieves the minimal theoretical sam- 
pling rate derived in [15]. In addition, while other works [20], 
[18], [19], [17], [22] provided reconstruction algorithms only 
for signals defined over a finite union of subspaces, here we 
provide a first systematic sampling and reconstruction method 
for signals in an infinite union of subspaces. 



V. Related Sampling Problems 

In the introduction, we outlined previous approaches to 
time-delay estimation. In this section, we explore in more 
detail the relationship between our sampling problem and pre- 
vious related setups treated in the sampling literature: sampling 
signals from a union of subspaces [15], [16], compressed 
sensing of analog signals [17], [18], [19], [20], and FR1 
sampling [23], [24]. 

A. Sampling signals from a union of subspaces 

A signal model which received growing interest recently is 
that of signals that lie in a union of subspaces [15], [16], [20], 
[18], [19], [22]. Under this model each signal x (t) can be 
described as [15] 

x(t)e\JS 7 , (62) 

7er 

where <S 7 are subspaces of a given Hilbert space and T is 
an index set. The signal x{t) lies in one of the subspaces 
S 1 , however it is not known in advance in which one. Thus, 
effectively, to determine x(t), we first need to find the active 
subspace, or the index 7. 

Our signal model, given by ([3), can be formulated as in 
(1621 . As described in Subsection IIII-AI once the time delays 
are fixed, each signal x (t) lies in a SI subspace spanned by 
K generators. Therefore, the set of all signals of the form (01 
constitute an infinite union of SI subspaces, where 7 is the set 
of delays r, which can take on any continuous value in the 
interval [0, T], and 5 7 is the corresponding SI subspace. 

In [15], [16] necessary and sufficient conditions are derived 
for a sampling operator to be invertible over a union of 
subspaces. For the case of a union of SI subspaces, [15] 
suggests a sampling scheme, similar to that used in [17] 
and in this paper, comprised of parallel sampling channels. 
Conditions on the sampling filters are then given in order 
to ensure reconstruction of the signals from its samples. In 
addition, the minimal number of sampling channels allowing 
perfect recovery of the signal from its samples is shown to be 
2K. This leads to a minimal sampling rate of 2K/T which 
is achieved by our scheme. However, in [15] no concrete 
reconstruction algorithms were given that can achieve this 
rate. Furthermore, although conditions for invertibility were 



B. Compressed sensing of analog signals 

The sampling problem in [17] also deals with signals that 
lie in a union of SI spaces and provides recovery algorithms. 
However in the setting of [17] there are finite number of 
possible subspaces, in contrast to our case, where there are 
an infinite number of possible subspaces. 

The signal model in [17] can be described in terms of N 
generating functions ae(t) as 

x{t)= E d * N at (* - nT) , (63) 

\i\=Kn£% 

where the notation \£\ — K means a sum over at most 
K elements. Thus, for each signal there are only K active 
generating functions out of N total possible functions, but 
we do not know in advance which generators are active. In 
principle, such signals can be sampled and recovered using 
the paradigm described in Section IHIl corresponding to N gen- 
erating functions. Indeed, any signal of the form d63l clearly 
also lies in the SI subspace spanned by the N generators (t), 
where some of the sequences di[n] are identically 0. However, 
this would require a sampling rate of N/T, obtained by 7Y 
sampling filters. Since only K of the generators are active, 
intuitively, we should be able to reduce the rate and still be 
able to recover the signal. The main contribution of [17] is a 
sampling scheme consisting of 2K filters that is sufficient in 
order to recover x(t) exactly. 

We can formulate our problem as a finite union of SI spaces 
of the form d63l if we assume that the K unknown delays are 
taken from a discrete grid containing N possible time delays. 
Under this assumption the generating functions in d63l can be 
expressed as 

ai (t) =g{t-t t ), 1 < I < N. (64) 

Therefore, under a discrete setting, the method of [17] can 
provide a sampling and reconstruction scheme for a signal of 
the form © with rate 2K/T. 

Similar to our approach here, the sampling scheme in [17] 
is based on 2K parallel channels, each comprised of a filter 
and a uniform sampler at rate 1/T. However, in order to 
achieve this minimum sampling rate, the reconstruction in 
[17] involves brute-force solving an optimization problem 
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with combinatorial complexity. The complexity of the recon- 
struction stage can be reduced by increasing the number of 
channels, which entails a price in terms of sampling rate. In 
contrast, our reconstruction algorithm is based on the ESPRIT 
algorithm, and can obtain the minimal sampling rate of 2K/T 
in polynomial complexity. Furthermore, we do not require 
discritization of the time delays but rather can accommodate 
any continuous set of delays. In this sense we can view our 
sampling paradigm as a special case of compressed sensing 
for an infinite union of SI spaces. Since previous work in this 
area has focused on sampling methods for finite unions, this 
appears to be a first systematic example of a sampling theory 
where the subspace is chosen over an infinite union. 

Another difference with the approach of [17] is the design 
of the sampling filters. In our method, we have seen that 
simple sampling filters can be used, such as low pass filter or 
bandpass filter-bank. In contrast, the scheme of [17] requires 
proper design of the sampling filters, which is obtained in 
two stages. In the first stage, N filters hi (t) , 1 < I < N 
are chosen that satisfy some conditions with respect to the N 
possible generating functions. At the second stage, a smaller 
set of p > 2K filters Si(t) is constructed from hi(t), via 

N N 

s *( t ) = Y,Y,Y, A ^ c *™ N hm (t ~ nT) , (65) 

l=\ 771=1 TlSZ 

where A is a p x N matrix that satisfies the requirements 
of compressed sensing in the appropriate dimension [36], and 
Ci m [n] are a set of sequences given explicitly in [17]. In order 
to arrive at filters that are easy to implement, a careful choice 
of the parameters is needed, which may be difficult to obtain. 

C. Signals with finite rate of innovation 

Another interesting class of signals that has been treated 
recently in the sampling literature are FRI signals [23], [24]. 
Such signals have a finite number of degrees of freedom per 
unit time, referred to as the rate of innovation. Examples of 
FRI signals include streams of diracs, nonuniform splines, and 
piecewise polynomials. A general form of an FRI signal is 
given by [23] 

aj(t) = 5^c n 0(t-t n ), (66) 

TiGZ 

where <j> (t) is a known function, t n are unknown time shifts 
and c„ are unknown weighing coefficients. Recovery of such 
signals from their samples is equivalent to the recovery of the 
delays t n and the weights c n . 

Our signal model (01 can be seen as a special case of 
(I66K where additional shift invariant structure is imposed. 
This means that in each period T the time delays are constant 
relative to the beginning of the period, whereas in a general 
FRI signal the time delay can vary from period to period. Our 
method is designed in such a way that it utilizes this extra 
structure to reduce the rate, while still guaranteeing perfect 
recovery. 

The FRI signals dealt with in [23], [24] are divided into 
three main classes: periodic, finite length and infinite length. If 
we address our signal model as an FRI signal it will generally 



fall into the third category of infinite length FRI signals. 
Some special classes of finite (and periodic) FRI signals where 
treated in [23], such as streams of diracs. For these special 
settings sampling theorems where derived with very specific 
kernels, that achieve the minimal rate (the rate of innovation). 
However, these methods are not adapted to the general model 
I®. 

Sampling and reconstruction of infinite length FRI signals 
was treated in [24]. The method in [24] is based on the use 
of specific sampling kernels which have finite time support: 
kernels that can reproduce polynomials or exponentials. In 
addition the function 4> (t) is limited to diracs, differentiated 
diracs, or short pulses with compact support and no DC 
component. The reconstruction algorithm proposed in [24] is 
local, namely it recovers the signal's parameters in each time 
interval separately. Naive use of this approach in our context 
has two main disadvantages. First, in our method the unknown 
delays are recovered from all the samples of the signal x (t). 
A local algorithm is less robust to noise and does not take 
the shared information into account. In addition, in terms of 
computational complexity, in our method all the samples are 
collected to form a finite size correlation matrix, and then the 
ESPRIT algorithm is applied once. Using the local algorithm 
requires applying the annihilating filter method, used for FRI 
recovery, on each time interval over again. 

A final disadvantage of the FRI approach is the higher 
sampling rate required. In order to discuss the sampling rate 
achieved by the local algorithm proposed in [24], we limit 
our discussion to the case where the function <fi (t) is a dirac, 
which is the main case dealt with in [24]. The theorems for 
unique recovery of the signal from its samples in [24] require 
that in each interval of size 2KLT S there are at most K diracs, 
where L is the support of the sampling kernel and T s is the 
sampling period. Since in each interval of size T we have K 
diracs, it can be easily shown that the minimal sampling rate 
is 2KL/T, which is a factor of L larger than the rate achieved 
by our scheme. For example, when using a B-spline kernel, 
which is the function with the shortest time support that can 
reproduce polynomials of a certain order, an order of at least 
N = 2K — 1 is needed, which has time support L = 2K. 
Thus, the sampling rate is 2K times larger than our approach. 



VI. Application 

In this section we describe a possible application of the 
proposed signal model and sampling scheme to the problem of 
channel estimation in wireless communication [37]. In such an 
application a transmitted communication signal passes through 
a multipath time-varying channel. The aim of the receiver is 
to estimate the channel's parameters from the samples of the 
received signal. 

We consider a baseband communication system operating in 
a multipath fading environment with pulse amplitude modula- 
tion (PAM). The data symbols are transmitted at a symbol 
rate of 1/T, modulated by a known pulse g (t). For this 
communication system the transmitted signal x t (t) is given 
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by 

N sym 

x T (t) = ^ d [n] g (t - nT) (67) 

n=l 

where d [m] are the data symbols taken from a finite alphabet, 
and N sym is the total number of transmitted symbols. 

The transmitted signal xt (t) passes through a baseband 
time-varying multipath channel whose impulse response is 
modeled as [38] 

K 

h(T,t)=J2°tk(t)S(T-T k ) (68) 

k=l 

where a k (t) is the path time varying complex gain for the fcth 
multipath propagation path and r k is the corresponding time 
delay. The total number of paths is denoted by K. We assume 
that the channel is slowly varying relative to the symbol rate, 
so that the path gains are considered to be constant over one 
symbol period: 

a k 0) = a k [nT] for t E [nT, (n + 1) T] . (69) 

In addition, we assume that the propagation delays are con- 
fined to one symbol, i.e r k € [0, T). Under these assumptions, 
the received signal at the receiver is given by 

K N sym 

X R ak N 9 (« - Th -nT)+n (i) (70) 

fc=l n=l 

where 

a k [n] = a k [nT] d [n] (71) 

and n (t) denotes the channel noise. 

The received signal xr (t) fits the signal model described in 
([3j. Therefore, if the pulse shape g (t) satisfies the condition 
( 1321 with p = 2K, our sampling scheme can recover the time 
delays of the propagation paths. In addition, if the transmitted 
symbols are known to the receiver, the time varying path gains 
can be recovered from the sequences a k [n]. 

As a result our sampling scheme can estimate the channel's 
parameters from samples of the output at a low rate, propor- 
tional to the number of paths. As an example, we can look 
at the channel estimation problem in code division multiple 
access (CDMA) communication. This problem was handled 
using subspace techniques in [39], [40]. In these works the 
sampling is done at the chip rate 1/T C or above, where T c 
is the chip duration given by T c = T/N and N is the 
spreading factor which is usually high (1023, for example, 
in GPS applications). In contrast, our sampling scheme can 
provide recovery of the channel's parameters at a sampling 
rate of 2K/T. For a channel with a small number of paths, 
this sampling rate can be significantly lower than the chip rate. 

Another example is UWB [41] communications which has 
gained popularity recently. In this technology the bandwidth 
of the transmitted pulse can be up to several gigahertz. 
Current technology commercial ADCs cannot operate at these 
sampling rates. For example, the highest sampling rate ADC 
device, manufactured by National Semiconductor, supports 
sampling rates of up to 3GHz at a relatively low resolution of 
8 bits and high power consumption. In contrast, our proposed 



method, has a potential of reducing the sampling rate, into 
rates which can be achieved by lower rate ADCs with better 
resolution and lower power consumption. 

VII. Numerical Experiments 

We now provide several experiments in which we exam- 
ine various aspects of our proposed method. The numerical 
experiments are divided into 4 parts: 

1) demonstration of a channel estimation application; 

2) evaluation of performance in the presence of noise; 

3) effects of approximation of the matrix using only 
a finite number of measurement vectors; 

4) effects of imperfect digital correction filtering, using 
finite length filters. 

In all the simulations, except for the one in IVII-DI we 
use the sampling scheme described in Section IIH-C.U which 
consists of a bank of ideal band-pass filters. We assume that 
the working band is T = [0, 3jrp] , and that the function g (t) 
has constant frequency response on that frequency range, i.e, 
G (oj) — T, uo G T . In order to improve the robustness to 
noise in the delay recovery stage, we use the total least-squares 
(TLS) version of the ESPRIT algorithm described in [1 1]. All 
the results are based on averaging 1000 experiments. 

A. Channel estimation 

In the first simulation we demonstrate a channel estimation 
application. We consider a time-varying channel of the form 
d68l l. with K = 4 paths. In order to simulate a time varying 
channel, the channel's gain coefficients a k [nT] are modeled 
according to the Jakes' model [42] as a zero-mean complex- 
valued Gaussian stationary process with the classical U-shape 
power spectral density. In such a model the varying rate of 
each gain coefficient depends on the maximal Doppler shift 
fd- In order to simulate a slow varying channel, relatively 
to the symbol rate 1/T, we used for each path a maximal 
Doppler shift of = 0.05/T. The energy of each time- 
varying path gain coefficient was normalized to (1/2) k+1 . 
The path delays were drawn uniformly in the range [0, T]. 
For the estimation N sym = 100 symbols were used where the 
data symbols are assumed to be known. The samples at the 
output of each of the sampling channels were corrupted by 
complex-valued Gaussian white noise with an SNR of 15dB. 

The number of sampling channels is taken to be p = 
5, which is only one more than the number of unknown 
delays. Although we have seen that 2K sampling channels are 
required for perfect recovery of every signal of the form (0), 
for some signals lowering the number of sampling channels 
is possible. Indeed, according to Proposition [2] for signals 
with dim (span (b [A])) = K, the minimal number of sampling 
channels required is K + 1. We will demonstrate that for this 
example, K+1 channels are sufficient. 

In Fig. [4] the original and estimated channels are shown. 
Since the gain coefficients of the channel are time-varying, 
only their averaged energy over time is shown in the figure. 
In Fig. |5j we plot the magnitude of the original and estimated 
gains of the first path versus time. From Figs. [4] and [5] it is 
evident that our method can provide a good estimate of the 
channel's parameters, even when the samples are noisy. 
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Fig. 4. Channel estimation with p = 5 sampling channels, and SNR=15dB. Fig. 6. MSE of the delays estimation versus SNR, for K = 2 and p = 4. 




observed MSE 
CRB 




15 20 25 

Channels number 



Fig. 5. Estimation of the time-varying gain coefficient of the first path, Fig. 7. MSE of the delays estimation versus the number of sampling channels 



p = 5, SNR=15dB. 



p, for K = 2 and SNR=10dB. 



B. Performance in the presence of noise 

In the next simulations we examine the effect of SNR and 
the number of sampling channels on the error in the delays 
estimation. We choose K = 2 close delays, <i = 0.4352T and 
< 2 = 0.521T. The sequences a k [n},k = 1,2, n= 1,2,. . . 100 
are finite length sequences with unit power chosen according 
to Jakes' model with f d = 0.05/T. 

Under the setting of the simulation, which consists of a 
pulse with constant frequency response and ideal band-pass 
filters, from d23l it can be verified that the sampling sequences 
satisfy the following relation in the time domain 



c [n] = N (r) b [n] , neZ, 



(72) 



The Cramer-Rao bound (CRB) for unbiased estimators of 
9k = ^f-tk from the data c [n], was derived in [43] for 
this data model. The TLS-ESPRIT algorithm, used for the 
delays estimation in our method, is known to be asymptotically 
unbiased [44]. Experimentally we verified that, under the 
simulation setup, the bias of the delays estimation is low 
enough for SNRs above 15dB. Therefore, in this range of 
SNRs, the CRB derived in [43] can give a lower bound on the 
MSE of the delays estimation (up to factor of 1^-), assuming 
our specific sampling scheme. 



In Fig. [6] the mean-squared error (MSE) of the time delays 
estimation is shown versus the SNR, when using p = 4 
sampling channels. For comparison we also plot the CRB. 
The figure demonstrates that our method achieves the CRB 
for SNR> 15dB, which is the range that delays estimation 
can be considered as unbiased. 

In Fig. [7] the MSE of the estimation of the time delays 
is shown versus the number of sampling channels, for a 
constant SNR of lOdB. The results demonstrate that the 
estimation error can be improved by increasing the number 
of channels. Therefore, oversampling improves the robustness 
of our method to noise. 

C. Effects of imperfect approximation of Udd 

Next, we investigate the influence of estimating the matrix 
Udd using only a finite number of measurement vectors d [n] . 
This number effects the total delay of our method, since 
reconstruction of the sequences et^ [n] is performed only after 
the unknown delays t are recovered. In Fig. [8] the MSE of the 
delays estimation is shown versus the number of measurement 
vectors used for estimation of Udd- A constant SNR of 
20dB and p — 4 sampling channels are used. Two cases 
are illustrated: in the first, the sequences ctfc [n] are taken 
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Fig. 8. MSE of the delays estimation versus the number of samples used, 
for K = 2, p = 4 and SNR=20dB. 
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Fig. 9. MSE of the delays estimation versus SNR for different lengths of 
digital correction filter bank approximations. 



according to the Jakes' model with parameter fd = 0.05/T 
and in the second case fd — 0.1 /T is used, which corresponds 
to sequences with faster variation rate. Fig. [8] demonstrates 
that the MSE depends on the variation rate of the sequences. 
Intuitively the faster the sequences vary, the more information 
each new measurement vector d [n] contains, improving the 
estimation of Rdd- In addition, it can be seen that using only 
50 measurement vectors, yields a reasonable estimation of the 
delays in the case of fd — 0.1/T. The same estimation error 
is achieved using 80 measurement vectors, when using slow 
varying sequences. This result can be further improved by 
increasing the SNR or the number of sampling channels. 

D. Effects of imperfect digital filtering correction 

In the next simulation we examine the effects of approxi- 



mating the digital correction filter bank W 



J«>T\ 



by finite 



length digital filters. The length of the filters affects the delay 
of our scheme. To demonstrate this point, we arbitrarily choose 
a sampling scheme composed of 3 non ideal band-pass filters 
with a frequency response given by 



St (w) 



(1 -0M) cos (u- ?ft) 



otherwise, 



where 



Tt = 



(73) 



(74) 



These filters satisfy the conditions of Proposition [T] and 
can model realistic sampling filters with non-flat frequency 
response. In this case a non trivial digital correction filter bank 
is required, whose coefficients are calculated using the inverse 
DTFT of W" 1 (e jujT ). 

In Fig. [9] the MSE of the delays estimation versus the SNR 
is plotted for different lengths of filters. At low SNRs the 
dominant error is caused by the noise, while for high SNRs the 
error is mostly a result of the correction filter approximation. It 
can be seen that a 49 taps filters provide a good approximation 
to the correction filter bank, resulting in a delay of 24 samples. 
When working at SNRs below 40dB, filters with 11 taps 
provide a reasonable approximation. 



VIII. Conclusion 

In this paper, we considered the problem of estimating 
the time delays and time varying coefficients of a multipath 
channel, from low-rate samples of the received signal. We 
showed that this problem can be formulated within the broader 
context of sampling theory, in which our goal is to recover an 
analog signal x(t) that lies in a SI subspace, spanned by K 
generators with unknown delays. This class of problems can 
be viewed as an infinite union of subspaces. 

We showed that if the channel has K multipath compo- 
nents, or equivalently, if the SI subspace is generated by K 
generators, than under appropriate conditions on the sampling 
filters, a sampling rate of 2K/T is necessary and sufficient 
to guarantee perfect recovery of any signal x(t). Here T is 
the transmission rate, or the period of the generators. We 
developed sufficient conditions on the generators and the 
sampling filters in order to guarantee perfect recovery at the 
minimal possible rate. To recover the unknown time delays, we 
showed that our problem can be formulated within the context 
of DOA estimation. Using this relationship, we proposed an 
ESPRIT-type algorithm to determine the unknown delays from 
the given low rate samples. Once the delays are properly 
identified, the time varying coefficients can be found by digital 
filtering. 

Besides the application to time delay estimation, the prob- 
lem we treated here can be seen as a first example of a system- 
atic sampling theory for analog signals defined over an infinite 
union of subspaces. Recently, there has been growing interest 
in sampling theorems for signals over a union of subspaces 
[15], [16], [20], [18], [19], [17], [22]. However, previous work 
addressing stability issues and concrete recovery algorithms 
have focused on finite unions. Here, we take a first step in 
the direction of extending these ideas to a broader setting that 
treats analog signals lying in an infinite union. 
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